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Abstract 

We review two recently developed efficient methods for calculating rate constants 
of processes dominated by rare events in high-dimensional complex systems. The 
first is transition interface sampling (TIS), based on the measurement of effective 
fluxes through hypersurfaces in phase space. TIS improves efficiency with respect 
to standard transition path sampling (TPS) rate constant techniques, because it 
allows a variable path length and is less sensitive to recrossings. The second method 
is the partial path version of TIS. Developed for diffusive processes, it exploits the 
loss of long time correlation. We discuss the relation between the new techniques 
and the standard reactive flux methods in detail. Path sampling algorithms can suf- 
fer from ergodicity problems, and we introduce several new techniques to alleviate 
these problems, notably path swapping, stochastic configurational bias Monte Carlo 
shooting moves and order-parameter free path sampling. In addition, we give algo- 
rithms to calculate other interesting properties from path ensembles besides rate 
constants, such as activation energies and reaction mechanisms. 
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1 Introduction 



Molecular simulation has become indispensable as a modern tool to gain in- 
sight in the kinetics of processes in complex environment by supplying detailed 
atomistic information that is not (easily) experimentally accessible. Using ei- 
ther classical or ab initio based atomistic force fields [1,2], techniques such 
as molecular dynamics (MD) [3,4] can model reactive events on a reasonable 
realistic level. In contrast to most experiments where kinetic properties such 
as the reaction rate are obtained by measuring the macroscopic population 
densities of reactant and product states over a long time (seconds), molec- 
ular dynamics simulations have to obtain good statistics with much smaller 
systems (usually ~ 100 to 100000 molecules) in the accessible time range of 
nanoseconds- microseconds using a time step of a few femtoseconds, as dictated 
by the molecular vibrations. This small timescale and system size limits the 
application to activated processes with relatively low barriers between reactant 
and product states. The computation of rate constants with straightforward 
MD becomes inefficient when the process of interest has to overcome a high 
activation barrier because the probability to observe a reactive event on this 
time- and system-scale decreases exponentially with the barrier height. The 
system will spend a long time in one of the stable states and occasionally 
jump -in relatively short time- to the other state. This separation of time 
scales results in two state kinetics: the exponential relaxation of the popula- 
tion densities [5]. 

The time-scale problem is traditionally solved by a two-step reactive flux 
method [6,7,8,9]. One first calculates the free energy as a function of a re- 
action coordinate describing the process. The transition state theory (TST) 
rate constant is then related to the probability to be at the maximum of the 
free energy barrier. This rate is only an approximation and the second part 
of the reactive flux methods computes the correction, the transmission coeffi- 
cient, by starting many fleeting trajectories from the top of the barrier [6,7,8,9]. 
However, the success of this method depends strongly on the choice of reaction 
coordinate. If the reaction coordinate fails to capture the molecular mechanism 
the corresponding transmission coefficient will be extremely low, making an 
accurate evaluation of the rate problematic if not impossible. For high dimen- 
sional complex systems, for instance chemical reactions in solution, or protein 
folding, a good reaction coordinate can be extremely difficult to find and usu- 
ally requires detailed a priori knowledge of the transition mechanism. Hence, 
TST based reactive flux methods will be ineffective for complex processes for 
which no prior knowledge is available. 

Chandler and collaborators [10,11,12,13,14] devised a method for which no re- 
action coordinate is needed, but only a definition of the reactant and product 
state. This method, called transition path sampling (TPS), gathers a collec- 
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tion of trajectories connecting the reactant to the product stable region by 
employing a Monte Carlo (MC) procedure called shooting and shifting. The 
resulting path ensemble gives an unbiased insight in the mechanism of the re- 
action. TPS has been successfully applied to such diverse systems as cluster 
isomerization, auto-dissociation of water, ion pair dissociation and on isomer- 
ization of a dipeptide, as well a reactions in aqueous solution (see Ref. [13] for 
an overview). A drawback of TPS is that the calculation of rate constants is 
rather computer time consuming. We therefore developed the more efficient 
transition interface sampling (TIS) method [15]. TIS allows a variable path 
length, thereby limiting the required MD time steps to the strict necessary 
minimum. The TIS rate equation is based on an effective positive flux for- 
malism and is less sensitive to recrossings. The shifting moves used in TPS 
to enhance statistics, are unnecessary in the TIS algorithm. Also, multidi- 
mensional or even discrete order parameters can easily be implemented in 
TIS. Recently, we showed that for diffusive processes one can exploit the loss 
of correlation along trajectories. This lead to the development of the partial 
path TIS (PPTIS) method, a variation of TIS that samples much shorter 
paths [16]. 

In this paper we re-derive the basic concepts of TIS and PPTIS in a more intu- 
itive way and relate them to the calculation of the transmission coefficient. For 
the mathematical validation of the expressions we refer to Refs. [15,16]. The 
paper is organized as follows. In section II we discuss the relation between 
several different microscopic expressions for the phenomenological rate con- 
stant present in the literature, and derive the positive effective flux formalism 
on which both interface path sampling methods are based. In section III we 
present the TIS and PPTIS formalism and precise algorithm. In section IV, we 
introduce new algorithms for alleviating ergodicity problems that might occur 
in path sampling simulations. The last section V, is reserved for new ways of 
extracting interesting properties from path ensembles, such as the activation 
energy of a reaction. We end with concluding remarks. 



2 Microscopic rate equations 

The calculation of reaction rate constants by computer simulation requires 
an expression for the rate constant in terms of microscopic properties. Such 
a microscopic rate expression needs a proper characterization of the reactant 
state A and product state B for each separate reaction, but should not be too 
sensitive to these state definitions, otherwise an unrealistic ill-defined rate will 
result. Once we have a rate expression, there are several ways to compute the 
reaction rate. The standard reactive flux method measures the flux through 
a single hypersurface in phase-space dividing the reactant state A from the 
product state B. In TPS the rate constant is taken from a time derivative of a 
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correlation function, which can be calculated by slowly confining a completely 
free path ensemble to an ensemble that connects reactant to product. The 
TIS approach measures a reactive flux through many interfaces between A 
and B. These three methods can be related to each other, as they ultimately 
compute the same properties. The TPS correlation function at t — becomes 
equivalent to the TST approach when A and B are adjacent in phase space 
[14] . The TIS effective positive flux formalism for a single interface is equivalent 
to TST-based transmission coefficient calculations [15]. The TIS rate equation 
can also be recast in terms of a TPS-like correlation function but then based 
on the so-called overall states of the system. In the following subsections, we 
will explain the reactive flux, TPS, and TIS methods and their connections in 
detail. 



2.1 Transition state theory 



The first step in TST is to choose a reaction coordinate A describing the 
transition from a stable reactant state A to a stable product state B. This 
reaction coordinate can be any function X(x) of phase space point x = {r, p}, 
with r the particle coordinates and p the momenta. Next, the free energy 
F(X) = — fcsTln(P(A)) is calculated by determining the probability -P(A) to 
be at A using, for instance, biased sampling techniques [17,18,19]. Here, ks is 
the Boltzmann constant and T is the temperature. The maximum A* in F(X) 
defines the dividing surface {x\X(x) = A*} separating state A from state B. 
By convention, the system is in A if X(x) < A* and in B if X(x) > A*. For a 
phase point x in A, the probability to be at the top of the barrier is: 



(S(X(x)-X*)) e-W> 
1 A, - fe - 1 - [9(X*-X(x)))- J^dXe-^y W 



where the brackets (. . .) denote the equilibrium ensemble averages, 6(x) and 
S(x) are the Heaviside step-function and the Dirac delta function, respectively, 
and P = (ksT)^ 1 . TST assumes that trajectories that cross A* do not recross 
the dividing surface Hence, the TST expression is equivalent to the positive 
flux through the dividing surface A*: 



k T A s B T = (X(x)6(X(x)))^P(X% eA , (2) 

where the dot denotes a time derivative and the subscript A* to the ensemble 
brackets indicates that the ensemble is constrained to the top of the barrier 
on the dividing surface A*. The TST rate constant is sensitive to the choice of 
reaction coordinate X(x) and will only be correct if the surface {x|A(x) = A*} 
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corresponds to the true transition state dividing surface: the so-called sepa- 
ratrix. For complex systems, it is impossible to know the location and shape 
of this curved multidimensional separatrix. It is possible, however, to correct 
the TST expression with a dynamical factor that is called the transmission 
coefficient. 



2.2 Transmission coefficients 



Traditionally, the dynamical corrected rate constant is derived by applying 
a small perturbation to the equilibrium state and invoking the fluctuation- 
dissipation theorem [20,5,4]. This leads, for instance, to the well known Bennett- 
Chandler (BC) [8,9] expression for the reaction rate 



= (\(x )6(\(x )-\*)e(\(x t )-\*)) 
XAB 'l(\*-\(x ))) 



where x t specifies the coordinates and momenta of the system at time t as 
obtained from a short molecular dynamics (MD) trajectory starting at Xq. The 
ensemble average is taken over all phase points x . For exponentially relaxing 
two state kinetics with a well defined rate constant, there is a separation of 
timescales: the reaction time r rxn (or expectation time for one single event) 
is much longer than the molecular time r mo i that the system spends on the 
barrier. In that case, Eq. (3) will reach a plateau value for r mo i <d t T rxn , 
which is equal to the correct phenomeno logical rate constant /cab- The function 
k^(t) will sensitively depend on the choice of the reaction coordinate A, but 
the plateau value will not. In the limit t — > + , the BC rate reduces to the 
TST expression Eq. (2) 

The transmission coefficient is defined as the ratio between the real rate con- 
stant and the TST expression: k = ^abI^/^b ^' 



K 



BC 



X(x )9(X(x t ) - X* 

(t) = ' ) n/' A * ( 4 ) 

X(x )9^X(x " 



The numerator in Eq. (4) counts trajectories with a positive but also with a 
negative weight. The latter trajectories leave the surface at t — with a nega- 
tive velocity X(xq), but are eventually found at the B side of the surface after 
a (few) recrossing(s). However, untrue B — > B trajectories do not contribute 
to the rate because the positive and negative terms cancel 1 (See Fig. 1). 
Similarly, the A — > B trajectories with multiple A* crossings are effectively 



This cancellation might seem to be not apparent if a trajectory recrosses the same 
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Fig. 1. Illustration of the difference in counting in the transmission coefficient 
Eqs. (4), (5), and (16). For simplicity, assume that the system consists of three 
kind of possible trajectories, as shown by this figure, that cross the dividing 
surface with the same speed v orthogonal to the surface. To the seven phase 
points on the surface (from top to bottom) the numerator of Eq. (4) will assign 
the values [— v, v, v, — v, v, 0, 0], while these are [0, 0, v, — v, v, 0, 0] for Eq. (5) and 
[0,0,^,0,0,0,0] for Eq. (16). The sum of these give the same result v. Evaluation 
of Eq. (16) in an actual computer algorithm requires the fewest MD steps as only 
phase points similar to the 3rd and 7th phase points would need the integration 
until reaching stable state regions. For instance, the fifth crossing point can be as- 
signed zero already as soon as one detects that its backward trajectory recrosses 
the surface. 

counted only once [15]. Although Eq. (3) gives the correct rate constant, it is 
rather unsatisfactory to sample only trajectories forward in time not knowing 
which contribute to the rate and which do not. Therefore, alternative expres- 
sions for the rate constant have been proposed taking the past into account. 
Here, they are referred to as the BC2 [8,9] expression 



;A(x„)e(A(x„))) 



and the positive flux PF [21] expression 



_ [x(x )e(x(x ))e(\(x t ) - a')) a . (\(x Q )e(\(x ))e(\(x_ t ) - x 



\{x Q )9(\{x Q ))) (\{x Q )6(\{x Q 



In Eq. (5) the theta functions guarantee that only true A — > B events are 
counted. Still, the numerator in Eq. (5) contains negative terms: those phase 
points Xq with a negative velocity X(x ) and with corresponding backward 



surface but with a different velocity. Still, this is the case. The absolute value of 
the flux of a trajectory is at each intersecting surface the same. A lower crossing 
velocity A is compensated by a higher probability to measure the crossing point as 
the trajectory spends more time at the surface. 



(6) 
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and forward trajectory that ends up in A and B, respectively. Eq. (6) counts 
only positive crossings, but cancellation with a negative term can occur when 
the backward trajectory also ends up at the B side of the dividing surface. 
At first sight, Eq. (6) seems to overcount A — > B trajectories with multiple 
A* crossings. However, if one realizes that each A — > B trajectory has an 
equivalent trajectory B — > A by reversing the time, an overall cancellation of 
positive and negative terms ensures a proper final outcome. 

For completeness, we mention that there are also similar expressions by Berne 
[22,23] and a relation by Hummer [24] that counts both positive and negative 
crossings with a positive weight, but only if the corresponding trajectory ends 
at opposite sides of the surface and with a weight lower than |A| if its trajectory 
has more than just one crossing. Ruiz-Montero et al. designed a transition zone 
method in which they measure the flux on many places at the top of the barrier 
and weight them to the inverse free energy [25]. 



2.3 The effective positive flux formalism 

A more intuitive, yet sound, alternative to the above expressions is the ef- 
fective flux formalism. We can illustrate this formalism with an analogy to 
the migration of people from country A to B. To determine the emigration 
rate we can simply count the number of persons that cross the border from 
A to B within a certain time interval. However, we should not count tourists. 
This group consist of people who have a nationality A and will only spend a 
short time in B, or have a nationality B and are actually on their way back. 
Moreover, we have to be aware that some emigrants might cross the frontier 
several times on their way. To prevent overcounting, we should only count one 
specified crossing for each person, for instance, the first or the last crossing 
of the emigration journeys from A to B. The same reasoning can be applied 
when calculating the rate constant of a reaction. In a molecular simulation 
we can check the 'nationality' of the system and the one-crossing condition 
by simply following the equations of motions backward and forward in time. 
This procedure, to count only true events and to avoid counting recrossings 
is what we call the effective positive flux formalism. In Sec. 2.5 we give the 
mathematical expression of the effective positive flux. 

It is surprising that the effective positive flux counting strategy is not so 
common. To our knowledge only two slightly different expressions of a trans- 
mission coefficient based on the effective positive flux have been proposed in 
Refs. [26,27]. In all other expressions found in the literature the counting of 
recrossings is not avoided, but the final rate constant follows through cancel- 
lation of many negative and positive terms. The effective flux transmission 
coefficients formulation is most useful when applying a single dividing surface 
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and when recrossings are apparent [23] . In general, we note that any averaging 
method counting only zero and positive values will show a faster convergence 
than one that is based on cancellation of positive en negative terms. Moreover, 
in the effective flux formalism many trajectories will be assigned as unreactive 
after just a few MD steps (See Fig. 1), thus reducing the number of required 
force evaluations. A comparative study of ion channel diffusion [23] showed 
that the algorithm based on effective positive flux expression of Anderson [26] 
was superior to the other transmission rate expressions. Moreover, it was found 
as efficient as an optimized version of the more complicated Ruiz-Montero 
method [25]. 



2.4 TPS correlation function 



In TPS one also has to define an order parameter X(x), but this does not 
have to be a properly chosen reaction coordinate capturing the essence of the 
dynamical mechanism. Instead, it is sufficient but necessary that this function 
is able to characterize the basins of attraction of the stable states [13]. By 
definition the system is in A if X(x) < Xa and in B if X(x) > X B with < 
Xb- Clearly, the two states are not connected and the intermediate barrier 
region, belongs neither to A nor to B. By introducing following characteristic 
functions 



Ha{x) = 1, if x E A, else Ha{x) = 

^b(^) = 1, if x e B, else h B {x) = 0. (7) 

the TPS-correlation function is defined as: 

r(A (h A (x )h B (x t )) 

C(t) ~ <M*o)> • (8) 

If there is a separation of timescales, this correlation function grows linearly in 
time, C{t) ~ kABt, for times r mo i < t < r rxn In that case, the time dependent 
reaction rate 

k T A P B S (t) = ±C(t) (9) 

reaches a plateau for r mo i < t < r rxn . C(t) can be calculated in a path sam- 
pling simulation employing the shooting and shifting Monte Carlo moves, in 
combination with an umbrella sampling algorithm in which the final region B 
is slowly shrunk from the entire phase space to the final stable state B [14]. 
The disadvantage of such a procedure is that it can take a relatively long time 
T moi before C(t) reaches a plateau (longer than in a transmission coefficient 
calculation [14]). 

All paths in the path sampling should have a minimal length T > r mo i and 
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as a result unnecessarily long periods are spent inside the stable state basins 
of attraction. Moreover, inspection of Eqs. (8) and (9) shows that a necessary 
cancellation of positive and negative terms can slow down the convergence of 
the MC sampling procedure. In the case of adjacent A and B regions, the TPS 
formalism becomes equivalent to the TST approximation in the limit t — > 
[14]. 



2.5 The road to TIS 



The TIS method is based on the measurement of the fluxes though multiple 
dividing surfaces. Consider a set of n + 1 non-intersecting multidimensional 
interfaces {0, 1 ... n} described by an order parameter X(x) that does not have 
to correspond to the real reaction coordinate. We choose Aj, % — . . . n such 
that Aj_i < Aj, and that the boundaries of state A and B are described by Ao 
and A n , respectively. For each phase point x and each interface i, we define a 
backward time t\{x) and forward time t{(x): 



t^xo) = - max [{t\X(x t ) = Aj A t < 0}] 

t{(x ) = + min [{t\X(x t ) = Aj A t > 0}] , (10) 

which mark the points of first crossing with interface i on a backward (forward) 
trajectory starting in x$. Note that t\ and t{ defined in this way always have 
positive values. Following Ref. [15], we then introduce two-fold characteristic 
functions that depend on two interfaces % ^ j, 



1 if tj(s) <<*(*), ^ |l if if (s) <«£(*) 

otherwise 1,3 [0 otherwise 

which measure whether the backward (forward) time evolution of x will reach 
interface % before j or not. However, as the interfaces do not intersect, the 
time evolution has to be evaluated only for those phase points x that are in 
between the two interfaces % and j. In case % < j, we know in advance that 
t b i' f (x) < tY(x) if X{x) < Aj and tf J \x) > tf (x) if X{x) > Xj. When the 
system is ergodic, both interfaces i and j will be crossed in finite time and 
thus h\j(x) + h) ti (x) = h{j(x) + hj ti (x) = 1. The two backward characteristic 
functions define the TIS overall states A and £>: 




M*) = K,n(x), h B (x) = h b nfi (x). (12) 

Together, the overall states cover the entire phase space and, within certain 
limits, do not sensitively depend on the precise boundaries of stable states A 
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and B. With these new characteristic functions we can write down a correlation 
function similar to Eq. (8): 



cm = <^2ftM> , (13) 

(h A (xo)) 



This correlation function exhibits a linear regime ~ k^t for < t < r r . 
Therefore, we can simply take the time derivative at t — yielding 



(h b (x )X(x )S(X(x )-X n ) 

kw= - ■ (14) 



One can easily verify that here only positive terms contribute to the rate. 
The connection to the transmission coefficient can be made by using following 
relation [15]: 



for Xi < Xj < Xk- Using this equality, we can write down a transmission 
coefficient similar to the ones in Sec. 2.2 but then based on the effective positive 
flux [27]: 



TTq (hl(xo)X(x )9(X(xo))hUxo) ; 

k = — — T . y. — l^- — ^ (i6) 



X(x )6(X(xo 



for Xi = X*. Although, in principle #^A(a;o)J is redundant in the numerator of 

Eq. (16) as h b Q ^(xq) = if X(xq) < 0, it is there to highlight that only positive 
crossings are counted. Trajectories started at x on interface % are followed 
backward in time until they reach stable region A or recross interface i. Then, 
only the ones that reach stable region A are also followed forward in time 
until they reach one of the stable regions. The slightly different effective flux 
expression of Ref. [26] follows trajectories until reaching the plateau region 
time and counts for each A — > B trajectory only the last crossing instead of 
the first. 
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3 (Partial path) transition interface sampling 

3. 1 Formalism 

In a system for which the correct reaction coordinate A is known in advance 
and that is not dominated by recrossings, the effective positive flux formal- 
ism of (Eq. (16) and Ref. [26]) is probably the best choice when using a single 
dividing surface [23] . However, for complex systems, for instance chemical reac- 
tions in solution, any intuitively chosen reaction coordinate can give arbitrary 
small transmission coefficients, making an accurate computation prohibitive. 
To improve reaction coordinates by e. g. taking solvent degrees into account is 
generally a difficult job. Some progress has been made by using the coordina- 
tion number as reaction coordinate [28,29], but this ad hoc approach probably 
only works for specific systems. For instance, we showed that a proton transfer 
reaction in water depends very sensitively on the angular orientation of the 
surrounding water molecules [30]. Similarly, the degrees of freedom in a protein 
are so large that dynamical folding processes are at best only very qualita- 
tively described by order parameters. Quantities such as radius of gyration or 
number of native contact do usually not correspond to reaction coordinates 
[31]. Subtle effects, e.g. the solvent structure, play also here a role. To incor- 
porate all these subtleties in a single one- dimensional reaction coordinate is 
an immense task and can only be successful if the precise reaction mechanism 
is already known in advance. The TPS and TIS techniques do not rely on a 
reaction coordinate. The TIS hypersurfaces do not have to coincide with the 
transition state dividing surface. 

At the end of this section we give TIS (and PPTIS) rate expressions that 
can be employed in a computer algorithm. First, as the derivation of the TIS 
and PPTIS formalism requires a proper notation, we introduce following flux 
function 

Mx) = h b jti {x)\X(x)\S(X{x) - Xi) = h b jti [x) hm -Le(At - t{(x)) (17) 

The first equality has the same flux notation as Eq. (14), but the second 
equality is more useful in practice. An MD trajectory might cross interface Aj, 
but consists of discrete time slices that are never exactly on the surface (as 
opposed to a transmission coefficient calculation). However, 4>ij{x) can still 
be defined for the discrete MD set of time-slices by taking At equal to the 
molecular time-step. In words, <pij{x) equals 1/At if the forward trajectory 
crosses Aj in one single At time-step and the backward trajectory crosses Xj 
before Aj. Otherwise (f>ij(x) vanishes. In addition, we introduce a flux function 
that incorporates also the forward trajectory 
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*$(x) = Mx)h{ tm (x) (18) 
By making use of Eq. (15) we can write for A, < Xj < A^: 



{Mx)) = (*$(x)) (19) 
and, thus, the rate constant (14) becomes 

k A B = (0n,O> / (h A ) = / (h A ) (20) 

for each Aj with < i < n. 

The second step is to define a conditional crossing probability that depends 
on the location of four interfaces: 



ptij) = (*$)/(<!>*)■ (2i) 

In words, this is the probability for the system to reach interface I before m 
under the condition that it crosses at t — interface i, while coming directly 
from interface j in the past (see Fig. 2). The probabilities in Eq. (21) are the 
building blocks for both TIS as PPTIS to construct expressions for the rate 
constant. The probabilities P( l m \j) are defined on any set of four interfaces. 
The case, where m = j = and m = j = n, is of special interest for TIS and 
will be annotated as follows 

TMA.M = P(i\i), VBfrAXi) = P(Un) (22) 

For PPTIS, two types of crossing probabilities are required: the one interface 
crossing probabilities 



X m ^i ^i 

Fig. 2. The conditional crossing probability P(J n |*) for a certain configuration of 
interfaces Aj, Aj,Aj, and A m . The condition |*-) is depicted by the arrow and the 
solid line for two phase points (the dots): from this phase point one should cross 
Aj in one single At time-step in the forward direction, and, besides, its backward 
trajectory should cross Xj before Aj. Two possible forward trajectories are given by 
the dashed line. The upper crosses X m before A;, the lower crosses A; as first. The 
fraction whose forward trajectories behave like the last case equals P( l m \))- 
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Pt — P (i-l\i-l)i Pf — P (i+l\i+l) 

pT = P(H\li), pi = P(tl\li), (23) 
and the long distance crossing probabilities 

P+^PQl), P-^Ptft 1 ). (24) 

Using these probabilities, the TIS rate constant can be written in terms that 
can be determined in a computer simulation [15] 



kAB = \ Pa(A w [Ai) 

n-l 

^(A„|A 1 ) = n^(A m |A i ) (25) 

The first factor is a flux and can be calculated by straightforward MD as 
Ai will be close to A (see Sec. 3.2). The second factor, the crossing probability 
^(AnlAi), is calculated using the factorization in Eq. (25) into probabilities 
Va(\+i\\) that are much higher than the overall crossing probability. These 
can be calculated using the shooting algorithm as will be explained in Sec. 3.3. 

For PPTIS the set of equations are as follows [16]: 



kAB = ——PZ, k B A = 77 v p n I 26 ) 



(h A ) n) (h B ) 

p- '' i- '- ' p- ~ 1 izl (07) 

^3 , 1= D- ' ^3 ~ I± 1 D- y l< ) 



Pf-l P J + -l p- _ Pj-l P 3 

Pf-l + Pj-l P j-l ' 3 Pj-1 + Pj-l P j-l 



The factor tt^c is identical to the TIS flux factor, whereas to obtain the re- 

verse rate ksA only a single extra factor is needed. The P+ and Pf are 

obtained via the recursive relations (27) once all single crossing probabilities 
of Eq. (23) are known. Starting with P± = Pf = 1, we can iteratively deter- 
mine (P, + , Pf) for j = 2, . . . until j = n. The one-hopping probabilities (23) 
can again be calculated using the shooting algorithm. The PPTIS formalism 
basically transforms the process of interest into a Markovian sequence of hop- 
ping events. Yet, if the dynamics is diffusive and the interfaces are sufficient 
far apart, the rate formalism (26) and (27) will be almost exact [16]. 
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3.2 The flux algorithm 



The flux factor is the effective flux through Ai of the trajectories coming 
from A (from A). This factor is most conveniently computed with the first 
two interfaces identical. Although ffi 1 ^ is not well defined for Ai = Ao, we can 
think that Ai = Ao + e in the limit e — > 0. In this way, the effective positive 
flux will be equal to the simple positive flux through Ai (trajectories cannot 
recross without re-entering A, hence, all crossings are counted.). Similarly, for 
the reverse rate ksA we can set A n _i = \ n — e. If Ai is chosen close enough to 
A the flux factor can be obtained by straightforward MD initialized in A and 
counting the positive crossings through Ai = A during the simulation run: 

<M = 1 N c (28) 

(h A ) AtN MD { } 

with At the MD time step, N MD the number of MD steps, and the number 
of counted positive crossings. To calculate the rate at constant temperature 
instead of constant energy, one can apply a Nose-Hoover [32,33,34,35] or An- 
dersen [36] thermostat. However, one should be aware that these thermostats 
do give the correct canonical distribution at a given temperature, but mod- 
ify the dynamics in an unphysical way. Hence, static averages (A(x)) will 
be correct, but time correlation functions (A(x )B(x t )) most likely not. As 
N+ ~ /fl(Ai - X(x ))9(\(x At ) - Ai)) is actually a correlation function over 
a very short time, this effect will be small. However, if necessary one can easily 
correct for this by explicitly counting only phase points x that in absence of 
the thermostat will cross Ai in one At time-step. Applying this correction is 
computationally cheap as it does not require any additional force calculations. 
In Appendix A we describe some possibilities for further optimization of the 
flux algorithm. 

3.3 The path sampling algorithm 

To calculate the conditional probabilities in TIS and PPTIS we use a path 
sampling algorithm [14]. However, there are some differences with the classic 
TPS algorithm. Most importantly, in (PP)TIS the path length is variable, 
which has a small implication for the acceptance criterion for the shooting 
move. In appendix B we derive this acceptance rule for arbitrary (stochastic 
or deterministic) dynamics. The main tools in the MC sampling of trajectory 
space are the shooting move and the time-reversal move [14]. In particular 
for PPTIS time-reversal moves can be quite effective. Shifting moves that 
enhanced statistics in TPS are not needed and even useless in (PP)TIS. 
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TIS algorithm: The quantity of interest in TIS is the crossing probability 
Pa(\+i\\) (or Pb(\-i\\) for the reverse rate constant ksA)- To calculate 
this probability by sampling in the \ interface ensemble one needs an initial 
path that starts in A (at A ), crosses the interface A; at least once, and finally 
ends by either crossing Ao or Aj+i. In general one can take simply a successful 
path from the previous Aj_i interface ensemble that reached Aj, and complete 
its evolution till reaching either A or Aj+i. For more details on initial path 
generation we refer to Ref. [14]). The phase space point x is then defined 
as the first crossing point of this path with interface Aj. It is convenient to 
use a discrete time index r = int(t/At), and let r h = \v&{t\{x^) j At) and 
= int(min[to(^o)i t{ + i(xo)]/At) be the backward and forward terminal time 
slice indices, respectively. Including xq, the initial path then consists of = 
r b + W + 1 time slices. Choosing a probability 7 < 1 and a Gaussian width 
a w we now start following loop: 

• Main loop 

(1) Take a uniform random number a,\ in the interval [0 : 1]. 

(2) If <y,\ < 7 perform a time-reversal move. Otherwise, perform a shooting 
move. 

(3) If the trial path generated by either the time-reversal or shooting move 
is a proper path in the Aj ensemble accept the move and replace the old 
path by the new one, otherwise keep the old path. Update averages and 
repeat from step 1. 

• Time-reversal move 

(1) If the current path ends at Aj+i reject the time-reversal move and return 
to the main loop. 

(2) If the current path starts and ends at Ao, reverse the momenta and the 
order of time-slices. On this reverse path, xq is the new first crossing point 
with Aj. Return to the main loop. 

• Shooting move 

(1) On the current path with length choose a random time slice r', with 

-T b < T' < T f . 

(2) Change all momenta of the particles at time-slice r by adding small ran- 
domized displacements Sp = 5w^/m with 5w taken from a Gaussian dis- 
tribution with width a w and m the mass of the particle [14]. 

(3) In case of constant temperature (NVT) simulations: accept the new mo- 
menta with a probability [4] : 



Here, E(x) is the total energy of the system at phase space point x. In case 
of constant energy (NVE) simulations in which possibly also total linear- 
or angular momentum should be conserved: rescale all the momenta of 
the system according to the procedure described in Ref. [37] and accept 
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the new rescaled momenta. 

If the new momenta are accepted continue with step 4, else reject the 
whole shooting move and return to the main loop. 

(4) Take a uniform random number a 2 in the interval [0 : 1] and determine a 
maximum allowed path length for the trial move by: 

N£L = mt(iV (o) /«2). 

(5) Integrate equations of motion backward in time by reversing the momenta 
at time slice r', until reaching either Ao, Aj+i or exceeding the maximum 
path length N^l x . If the backward trajectory did not reach Ao reject and 
go back the main loop. Otherwise continue with step 6. 

(6) Integrate from time slice r' forward until reaching either A , A i+ i or ex- 
ceeding the maximum path length N^.. Reject and go back to the main 
loop if the maximum path length is exceeded or if the entire trial path 
has no crossing with interface \. Otherwise continue with the next step. 

(7) Accept the new path, reassign x to be the first crossing point with Aj and 
return to the main loop. 

Finally, the probability Pa(K+i\K) follows from: 

iMAm|Ai) " iytotai) (29) 

with N p (0 — > i + 1) the number of sampled paths that connect Ao with Aj+i 
and N p (total) the total number sampled paths in the ensemble of interface Aj. 

Time reversal moves do not require any force calculations. On the other hand 
two subsequent time reversals will just result in the same path. Therefore, 
we usually take 7 = 0.5 giving shooting and time reversal move an equal 
probability. Similar reasoning is applied to the choice of a w . If a w is large, 
many trial moves will fail to create a proper path. On the other hand a too 
small value of a w will result in almost the same path. Practice has shown that 
an optimal value of a w is established when approximately 40% of the paths is 
accepted [12]. This will usually imply that a w will be larger for the interfaces 
Aj close to A than the ones closer to B. The mass weighted momenta change 
at step 2 of the shooting algorithm is chosen such that the velocity rescaling 
at step 3 maintains detailed balance [37]. In principle, NVT simulations do 
not require rescaling and 5p can be taken from any symmetric distribution. 
The integration of the equations of motion at step 5 and 6 of the shooting 
move are normally performed by constant energy MD simulations without 
using a thermostat to describe the actual dynamics as realistic as possible. 
The temperature only appears at the acceptance criterion at step 3. In this 
algorithm we go from one phase point Xq to a new one by means of 
many MD steps. Therefore, it has a strong similarity with hybrid MC [38]. 
Hence, the argument that the dynamics should be time reversible and area 
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preserving [4] should also be applied here. For this reason, we strongly advice to 
use the velocity Verlet [39] algorithm rather than higher order schemes such as 
Runga-Kutta. The maximum allowed path length N^. in step 4 is introduced 
to maintain detailed balance when sampling paths of different length and to 
avoid having to reject very long trial paths afterward [15]. 

PPTIS algorithm: the four one-interface probabilities pf,p^,pf, and p\ for 
a single interface \ can be calculated simultaneously [16] with paths that start 
at Aj_i or Aj + i and end by crossing either Aj_i or Aj+i. All paths should have 
at least one crossing with \. Hence, r h = int(min[^_ 1 (xo), t^ +1 (x )] / At) and 
= int(min[tf_ 1 (x ), t{ +1 (x )] / At) . The path sampling is then identical to 
the TIS algorithm except that Aj_i is used instead of Ao, time reversal moves 
are always accepted and the backward integrating at step 5 is not rejected 
when reaching X i+1 as paths may start from both sides. The one-interface 
crossing probabilities are then given by 



N p (i - 1 



i + l) 



N p (i - 1 -> % + 1) + N p {i — 1 — > i — 1) 



Pi 



N p (i + l^i- 1) 



N p (i + 1 -> % - 1) + N p (i + 1 -> % + 1) 



pf, 



p i = i - pi 



(30) 



3.4 Defining the interfaces 



The order parameter A in TPS and TIS does not have to correspond to a 
reaction coordinate that captures the essence of the reaction mechanism. The 
only requirement is that A can distinguish between the two basins of attraction. 
In TIS this occurs via the two outer interfaces A and A n that define state A and 
B. The definitions of A and B are more strict than in TPS [15]. The boundaries 
Ao and A n should be defined such that each trajectory between the stable states 
is a rare event for the reaction we are interested in. In addition, the probability 
that after this event the reverse reaction occurs shortly thereafter must be as 
unlikely as an entirely new event. In other words, a trajectory that starts in A 
and ends in B is allowed to leave region B shortly thereafter, but the chance 
that it re-enters region A in a short time must be highly unlikely. Sometimes 
it is not sufficient for a proper definition of the boundaries A and A n to only 
use configuration space. In the dimer study of Ref. [15] an additional kinetic 
energy constraint was introduced to ensure the stability of state A and B. 

The intermediate interfaces can be chosen freely and should be placed to 
optimize the efficiency. This is, of course, system dependent, but reasonable 
estimates can be made a priori. Let us write down the total computation 
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time as CPU ~ Z)i=a NiLi with Nw the number of windows (interface en- 
sembles), Ni the number of paths in the ensemble of interface % required to 
obtain a desired precision e^, and Lj the average path length. Here, we ne- 
glect the influence of rejections and the fact that two successive pathways 
in the MC sequence are not completely uncorrelated. We chose the inter- 
face separations and the number of paths such that P(Aj+i|Aj) = p and 
Ni = n p , resulting in q = e for all %. The total error e tot , that we fix, is 
related by e^ ot = iYVe 2 with e 2 ~ (1 —p)/{pn p ). Hence, n p ~ Nw(l —p)/p 
yielding CPU ~ J2i=[ L>iNw(l ~ P)/P- The number of windows follows from 
pN w _ p(\ n \\ Q } =v. ]\[ w r-j — l/ln(p). Except for diffusive barrier cross- 
ings [16], that are most conveniently treated by PPTIS, the average path 
length Li has a linear dependence ~ i(X n — \q)/Nw [15]. Taking this all into 
account, the final result gives CPU ~ ln(p)~ 2 (l —p)/p that has a minimum for 
p = 0.2. Although, we made several assumptions in this derivation, we believe 
that in general P(A i+ i|Aj) « 0.2 for all % is close to an optimum efficiency 

Between the interface positions one can use of a finer grid of sub-interfaces to 
obtain the crossing probability function Pa( A|Ai) [15] which is the path space 
analogy to a Landau free energy profile P(A). For PPTIS different require- 
ments exist for the position of interfaces. As the PPTIS formalism is based on 
a complete memory loss over distances larger than the interface separations, 
the PPTIS interfaces should be set sufficiently far apart. The calculation of 
memory loss functions can help to determine the minimum required distance 
to establish this [16]. 

We would like to stress that although PPTIS transforms the system into a 
(pseudo) Markovian hopping sequence based on local transition probabilities, 
it still maintains considerable history dependence. For example, the chance 
to go from interface i to interface i + 1 is assumed to be equal for the path 
that arrived at i via the sequence i — 2 ^ i — 1 ^ i or via the sequence 
% — > % — 1 — > %. However, this transition to % + 1 from % can still be different 
when its history had hopping sequence i + 1 — > i. 



4 Improving the sampling 

4-1 Parallel path swapping 

Biased sampling methods such as constrained dynamics [19], multicanoni- 
cal [40] or umbrella sampling [17,18] can suffer from substantial ergodicity 
problems when the order parameters are not equal to the reaction coordinate. 
This lack of ergodicity usually shows up in hysteresis in the free energy curves 
(see e. g. [30,41]), and gives, besides a low transmission coefficient, rise to an 
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Fig. 3. Path swapping move for PPTIS. The last half of the path in the Aj ensemble 
and the first half of the path in the Aj+i are swapped to the Aj+i and Aj ensembles, 
respectively. 



additional error in the rate constant estimate. 



Transition path sampling was precisely devised to avoid this problem with 
reaction coordinates, and, in a way, also avoids ergodicity problems due to the 
non-local nature of the shooting move. This advantage showed up in the water 
trimer study [37] where the TPS algorithm was capable of finding two reaction 
mechanisms across different saddle points separated by a barrier higher than 
the total energy of the NVE simulation. We stress that this would have been 
much more difficult to achieve or even impossible in an umbrella sampling al- 
gorithm with several narrow windows. However, path sampling can also suffer 
from ergodicity problems if large barriers separate multiple reaction channels 
in a high dimensional rough energy landscape. In particular in the case of 
PPTIS, the short paths are much less likely to overcome such barriers. 

Parallel tempering techniques (also known as Replica Exchange methods) can 
facilitate the sampling [42], but requires a rather large computational effort 
and cannot be applied at constant energy. Here, we propose a less expen- 
sive parallel method especially tailored for PPTIS to enhance ergodicity. This 
parallel path swapping (PPS) technique is based on the exchange of paths 
between two subsequent interface ensembles. Fig. 3 shows one path in the Aj 
ensemble, consisting of all possible paths crossing Aj while starting and end- 
ing at either Aj_i or Aj+i, and one in the Aj+i ensemble consisting of all paths 
crossing Aj+i at least once, while starting and ending at either Aj or Aj + 2- We 
introduce a new MC move that attempts swapping the current path of the 
Aj ensemble with that of the A i+ i-ensemble, as depicted in Fig. 3. The swap 
move will be rejected if the Aj ensemble path does not end at or if the 
Aj+i ensemble path does not start at Aj. Otherwise, the move is accepted and 
the two trajectories are swapped from one ensemble to the other. Integrating 
the equations of motion backward (for the Aj ensemble) and forward (for the 
Aj +1 ensemble) will result in two entirely new paths for both ensembles. The 
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acceptance/rejection criterion appears before any expensive computation of 
MD trajectories. Moreover, once accepted we obtain a new path for both en- 
sembles for price of effectively only one path. This makes the path swapping 
move useful even if for systems not suffering from ergodicity problems. 

Another advantage of PPS is that it allows to go beyond the pseudo-Markovian 
description of PPTIS. Fig. 3 shows that the paths at the right hand side, if we 
include the dashed trajectory part, can connect four interfaces instead of only 
three. This extension allows for a long range verification of the memory loss 
assumption. Also, the development of new, smart algorithms based on PPS 
might be able to correct for memory effects or to search for ideal interface 
positions on the fly. 

While PPS is very effective when the confinement of short paths in PPTIS can 
cause sampling problems, even TIS and TPS algorithms might benefit from 
path swapping when multiple reaction channels exist. 



4-2 CBMC based shooting moves 

Originally developed to sample polymers at high densities, the Configurational 
Bias Monte Carlo (CBMC) technique grows chain molecules in a biased fashion 
in order to avoid unfavorable overlap of the beads [43,44,45,46]. The similarity 
between growing polymers and generating dynamical trajectories was the in- 
spiration for the development of TPS and has been exploited in the sampling 
of the stochastic path action [10,47]. However, this CBMC-like technique was 
found to be less effective than the shooting algorithm [14]. Here, we propose 
a combination of the shooting move with CBMC for diffusive systems that 
suffer from low acceptance due to a non flat rough free energy barrier. When 
shooting from one basin of attraction in such systems, the Lyapunov insta- 
bility causes the paths to diverge and return to the same basin of attraction 
before crossing the barrier. The use of some stochastic noise allows shooting 
in only one time direction and alleviates this problem slightly [48,31], but at 
the price that independent pathways are generated only after a number of ac- 
cepted shooting moves from the barrier region. This slow exploration of path 
space is even worse for processes proceeding via multiple dynamical bottle- 
necks, for instance reactions taking place though a short lived intermediate 
state 

Within the shooting algorithm, CBMC can be applied both at the shoot- 
ing point (the random time slice for which we change the momenta) and 
along the path by introducing some stochastic noise. At the shooting point 
t' we generate a set of N s momenta displacements {Sp^}, and accept these 
displacements using step 3 in Sec. 3.3. Each phase point is then integrated 
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Fig. 4. CBMC shooting move. (A) At the shooting point at the old trajectory 
(dashed line) four trial segments are released. In this example the momenta of 
segment 1 and 2 have been rejected and are not integrated further. Segment is 
retracing the old path. Of the trial segments, segment 3 has come farthest in its 
forward (solid arrow) and backward (open arrow) time evolution and will conse- 
quently have the highest weight. (B) The use of stochasticity allows the creation 
of trajectory jets at several points along the path. At each junction the path will 
follow the most favorable direction (bold solid line). The creation of trajectory jets 
at the old path (bold dashed line) is required to maintain super-detailed balance. 

forward and backward for a time tl, resulting in N s trajectory segments 
Sj = {x { fi_ TL)At , x$ +TL)At }, for j = 1, ... , N s (See Fig. 4-A). The time in- 
terval tl should be large enough to decide whether a trajectory has a chance of 
being successful, but much smaller than the average path length of a complete 
trajectory. All path segments are given a weight Wj 



w f = y(5p^)F( Sj ) (31) 

where \I> equals 1 (else 0) for accepted momenta changes Sp at the shooting 
point t'. The biasing function JF should be chosen to give the highest weight 
Wj to those segments that are most likely to produce a complete path of the 
corresponding interface ensemble. One possibility is to choose T = exp(aAA) 
with AA = \(x( T '+T L )At) ~ ^( x (r'-T L )At) and a a parameter optimized to the 
steepness of the barrier at x T i At . In that case, T is a function only of the 
backward and forward end points of the path segments Sj. The Rosenbluth 
factor for the set of trajectory segments is 



W*0 = 5>i (32) 
j'=i 

One of the segments Sj is selected with a probability Wi/W^ n \ To correct for 
this bias and to obey detailed balance, we also have to calculate the Rosenbluth 
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factor for the old path. The procedure is the same as above, but now we 
apply N s — 1 new random momenta changes {5p^} to the momenta of Sj at 
the same shooting point and again generate a set of segments of length 2tl At. 
This set is completed by adding segment s of the same length from the old 
path. The Rosenbluth factor for the old path equals 

N s -1 

W {o) = J2 w f- ( 33 ) 

3=0 

where is the weight of segment s , and with j = 1, . . . , N s — 1 are 
the weights for the segments that follow from {Sp^}. 

By imposing super detailed balance [4] the acceptance probability of segment 
% becomes 



-P aC c(so -> Sj) = min 



(34) 



Here, the weight functions Wi and the distributions p are still present, because 
they do not cancel as in the standard CBMC expression. The accepted segment 
is integrated to the complete path just as in the normal shooting move of 
Sec. 3.3. Of course, this procedure is computationally more expensive than the 
standard shooting move. However, the biasing function T allows to choose a 
segment with much higher probability to become a accepted path. We expect 
an increase in sampling efficiency when the gain in acceptance outweighs the 
cost of the construction of the trajectory segment sets. 

In the above algorithm we only can bias the growth of the first segment of the 
trajectory (the analog of the polymer in standard CBMC) because the rest of 
the trajectory follows deterministically once the first segment has been chosen. 
In the standard polymer CBMC a bias is introduced at each segment, and we 
can make use of the full power of CBMC if we consider stochastic trajectories. 
Introducing a small amount of stochasticity by for instance the Andersen 
thermostat [48] or by making use of the periodic boundary condition [49] will 
hardly change the dynamical properties of the transition process. 

Stochasticity allows us to create trajectory jets at several points along the 
paths (See Fig. 4-B). The first segment is created as in the deterministic 
procedure above. However, the chosen segment is not integrated to the full 
path length. Instead, we start with the end point of the forward trajectory 
and integrate a 'jet' of forward trajectory segments each evolving differently 
according to its own random noise. Each segment j of this 'jet' k has a weight 
Wjk similar to Eq. (31) and each jet will have a total weight Wk = J2j w jk- We 
select a segment i according to its relative weight w^/Wk, and continue with 
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the next jet of forward segments. The same is done for the backward paths, 
until the path is completed. After generating the new path, we have to repeat 
the 'jet' procedure for the old path as depicted in Fig. 4 in order to calculate 
the Rosenbluth factor of the old path. The total Rosenbluth factors are now 



^ = IW J . (35) 
k k 

where k runs over all the jets including the one at the shooting point avAt- 
The final acceptance criterion obeying super detailed balance is then 



P acc {o -> n) = min 
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(36) 



where is the segment weight at jet k on the old path and w^, is the 
weight of the selected segment of jet k on the new path. To take into account 
the change in path-length one should include a factor min 

[1, jV(°)/7V(n)] ; but 

this is usually implemented by defining a maximum path length as explained 
at step 4 of the shooting algorithm in Sec. 3.3. The above algorithm could be 
useful when the standard shooting move suffers from extreme low acceptance 
ratios. 



4-3 Time as transition parameter 



In TIS the choice of the order parameter is not critical as A does not have to 
correspond to the reaction coordinate. Yet, it is possible that the order pa- 
rameter A can bias the outcome of transition mechanism and rate constants - 
although much less than for the TST reactive flux method-, for instance, when 
the reaction mechanism leads in a direction that A does not allow. In principle, 
an order parameter-free sampling method is, therefore, highly desirable when 
examining unexpected contra-intuitive reaction mechanisms. One possibility 
for such a bias-free method is by using the time on the path outside A as 
transition parameter (we use 'transition' instead of 'order' to indicate that it 
is not a traditional order parameter as it is not based on a phase point). 

For a particular stable state A definition A , Pa(%+i\%) is the probability 
that a path, starting from A and remaining outside A over a time %, remains 
even longer outside A until at least T i+ i > %. To calculate the probability 
Pa{%+i\%) by a bias-free TIS simulation we generate an ensemble of trajecto- 
ries that have path lengths between % and %+i using the shooting algorithm 
of Sec. 3.3. At the shooting point, we integrate backward until reaching Ao or 
until the length of the trial trajectory exceeds T i+ i (or N^^At as defined at 
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step 4 of the shooting algorithm). If the backward trajectory exceeds either 
T i+l or N^ x At the shooting move is rejected. The forward trajectory is con- 
tinued until reaching A , or until a path length of T i+1 , or N^^At. The trial 
path is rejected if N^ x At is exceeded or if the trajectory ends at A in a time 
shorter than %. In the subsequent ensemble, the probability Pa (^+2 for 
T i+2 > %+i is calculated for all paths with at least a length 

This method, as illustrated in Fig. 5, will thus explore automatically the re- 




Fig. 5. Time as transition parameter. The square denotes the definition of the bound- 
ary for state A. The thin lines are free energy contour lines. The four panels show 
the representation of generated trajectories in successive time-interface ensembles. 
At panel 1), Pa{%+i\%) is the fraction of of trajectories that stay outside A longer 
than Ti+i (open arrows). All trajectories have at least a length %. The solid arrows 
are the paths that return to A before %+\. At panel 2), PaC^+2|^+i) is calculated 
for paths that remain outside A longer than % + \. The minimum length of the paths 
is further increased at panel 3). Incidentally, a path will end up in the yet unknown 
state B. At panel 4) the minimum path length constraint forces all the paths into 
the metastable state region B. From here, they will not return. Hence Pa(T\0) will 
show a plateau. 

gions further and further outside A. At some moment it will find the closest 
stable state region (state B). Trajectories reaching this region will not go 
back to A, hence, the overall crossing probability function P A (T\0) will show 
a plateau at some time T similar to standard TIS. 

Two-ended path sampling methods, such as TIS, PPTIS and TPS can only 
treat processes in which both stable states A and B are known. They cannot 
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find the final state starting from a single stable state, a fact already discussed 
by Dellago and Chandler [50]. The algorithm described here, might be a solu- 
tion to this problem. 



5 Extracting information from path ensembles 

5.1 Reaction mechanism 

The ensemble of paths collected by the TIS algorithm can be used to investi- 
gate the reaction mechanism. We believe that for this purpose the TIS path 
ensembles might even be more useful than the TPS path ensembles. The TPS 
method, first samples paths that all successfully reach B in the part to obtain 
the reactive flux function and then in the second step samples artificially short 
trajectories of fixed length to calculate the time correlation function C(t). Be- 
cause of this constraint, the resulting ensembles do not give useful information 
about the reaction. The TIS Aj-ensembles, on the other hand, contain the cor- 
rect distribution of paths that have crossed \ and are either going on to Aj+i 
or return to A. Some hidden order parameters can only be discovered by care- 
fully comparing configurations along reactive and unreactive trajectories that 
are similar in terms of order parameters which at first sight were considered as 
being the (only) important ones. For instance, the comparison of reactive and 
unreactive geometries with an almost identical orientation of the reactants 
showed that precise tetrahedral ordering of the solvent water molecules was 
an important factor in the hydration reaction of ketones [30]. Although there 
is currently no systematic way to extract the reaction coordinates from a path 
ensemble, once a reaction coordinate is postulated based on physical insight 
it can be tested using committor distributions [13]. 

5.2 Activation energies 

The activation energy E a is an important experimentally accessible quantity 
and is defined by the Arrhenius law 

k = Ae- pE % (37) 

where A is a system dependent prefactor. In fact, A and E a may also be 
temperature dependent. Such non Arrhenius behavior can be quite severe: 
sometimes reaction rates are even decreasing with increasing temperature, 
resulting in a 'negative activation energy' (see e.g. [51]). ^From Eq. (37) it 
follows that 
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An algorithm to calculate E a in a TPS simulation was given in Ref. [52]. 
Here, we use a similar approach to calculate E a in a canonical TIS simulation. 
Substitution of Eq. (25) in Eq. (38) results in 
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For any function A(x) we can write 
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with = (A(x)S(a;)) / (A(x)). Using 
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(E(x))^o = (E(x))^ +i o , 
most terms in Eq. (39) cancel, only leaving 
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E a = (E(x))^o - (E(x)) hA 

n — 1 , U - n 



(42) 



which is the difference between the average energy of state A and the energy 
of the transition pathways connecting A with B. Consequently, the calculation 
of the E a does not require all interface ensembles, but only the last ensemble 
A n _i. However, if all the path ensembles i = 1, ... ,n — 1 are available an 
activation energy function 



E a {\) = (E(x))^o i o - (E(x)) hA (43) 

can be calculated that should converge to a plateau analogous to the crossing 
probability P(A|Ai). A finer grid of sub-interfaces can be applied to obtain a 
continuous smooth function E a (\). 

Again, there is a subtle difference between the TPS and TIS algorithms. For 
the reaction rate determination, TPS requires a plateau in the time correlation 
function of Eq. (8), while TIS should give a plateau in A for the crossing 
probability P(A|Ai). Similarly, the TPS activation energy is expressed as a 
time dependent function that will converge to a plateau at times t = T [52], 
while the TIS activation energy reaches a plateau in terms of A. 
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Not all reactions show Arrhenius behavior. Therefore, it would be interesting 
to determine /cab(/3) for a range of temperatures. One can estimate the rate for 
a slightly different temperature by reweighting the crossing probabilities[17,18]. 
If (A(x)) is the average of an observable A(x) at inverse temperature (3, 
then (^A(x)e~ Af3E( - x ^ / ^e~ A/3B(x ^ should be the average at inverse temper- 
ature p + A(3. This reweighting technique can also be applied to the crossing 
probabilities (21) and the flux (28). Of course, Aj3 should be small to obtain 
good enough statistics. The calculation of the temperature dependence of in- 
dividual crossing probabilities has the advantage that the origin of possible 
non-Arrhenius behavior might be located (in terms of A) along the reaction 
path. 



6 Summary and Conclusions 



We reviewed the basic concepts of TIS and PPTIS and explained their rela- 
tion to TST based methods and TPS. We believe that path sampling meth- 
ods, TPS, TIS and PPTIS, are powerful when dealing with high dimensional 
complex process for which is a reaction coordinate is lacking. Among these 
methods, TIS can be considered as an improvement upon the original TPS 
giving a complete non-Markovian description of the reactive event, but more 
efficient. PPTIS improves the efficiency even more, but relies on the assump- 
tion of memory loss between interfaces. Hence, it should only be applied for 
diffusive barrier crossings. In addition to this review, we have introduced sev- 
eral new techniques in this paper. These novel methods comprise the CBMC 
based shooting moves, order parameter free methods, parallel path swapping 
and the calculation of activation energies. The efficiency of these methods 
should be tested by future simulations. We plan study this in the near future. 
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A The flux revisited 



In some cases, we can improve the efficiency of the flux calculation by sep- 
arating the flux into the probability to be on the Ai surface times a factor 
integrating over all possible velocities when leaving the surface. The flux term 
can then be calculated by combining straightforward MD with, as soon as we 
cross the Ai surface, the sampling of sets of randomized Gaussian distributed 
velocities A, after which the MD trajectory is continued with the old original 
momenta. In this way, we make optimal use of the statistics of the crossing 
points. The velocity sampling does not require force calculations and is there- 
fore cheap. In the following we assume that we always take Ai = Ao- Similar 
to Eq. (2) we can write: 



0(A)A\ P^Ua (A.l) 



(01,O) 

(h A ) \" V ' V "/Ai 

with 



p (AiW s ttM (A . 2) 

The two terms -^y- and P(\i) xe A can be obtained in the same MD simulation. 
As (S(X(x) — Ai)) dX is equal to the probability to find the system in the 
interval [Ai — \dX : Ai + \d\\, it can be measured by defining a width dX and 
performing a MD (or MC) simulation starting in A: 



1 N\ 

P(X 1 ) xeA = tt4t^ ( a -3) 
v UxeA dXN MB v ; 

with N\ 1 the number of counts in the specified interval and Amd the number 
of MD steps. However, this number can depend sensitively on the choice of 
bin width dX. Ideally one would like dX to be as small as possible, at the cost 
of having to perform a very long simulation run for a statistically accurate 
number N\ 1 . A better option is to weigh the crossings with a function de- 
pending on the velocity. Assume that we cross Ai in one MD step from xn± t to 
If dX is small neither of these points will lie inside the interval. How- 
ever, assuming a linear dynamics between these points, the system traverses 
from Xi&t to X(j + i)At in A su b equidistant sub steps. The number of phase points 
N\ 1 that lie in the dX interval of this short linear trajectory is approximately 
^AA sub /|A(a;(j +1 )Ai) — A(xjAt)|- The total number of MD moves A MD , of course, 
also increases by a factor A sub . So Eq.( A. 3) becomes 
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N M Di |A(X( i+ i)At) - HXiAt)\ 

where the * indicates that the summation has to be performed only for points 
i along the trajectory for which x iAt ^(i+i)At showed a crossing (positive or 
negative) with interface Ai. Further optimization can be achieved by writing 
|A(x (i+ i )At ) - A(>, At )| = |A x . At |A* + 0(At 2 ) = \X X(i+1) JAt + 0(At 2 ), but the 

velocity A at the interface would give the most exact result. If we also assume 
a linear change in time for the velocities between i and i + our best estimate 
for PiX^xeA is: 



p(\iUa = — ^S* tt7 1 x M ( A - 5 ) 

A MD Ai i \X{x iAt ; Ai)| 

where we have introduced the notation g(x iAt ; Xj) to denote the function g(x) 
at the crossing point of interface Xj obtained by a linear interpolation of the 
function between two successive trajectory points x iAt — > X(i + i) At : 



g(xiAt, Aj) = ^ — j _ ^ ^ j {[A(x^ + i)At) - \j]g(x iAt ) + [Xj - X(x iAt )]g(x (i+1)At ) 

(A.6 



The factor (9(X)Xj^ in Eq. (A.l) can be calculated in the same MD simu- 
lation with an additional sampling procedure. In some simple cases, there is 
even an analytically expression. For instance, in case the x-coordinate of par- 
ticle j is the order parameter, X(x) = Tj X in a constant temperature (NVT) 
simulation, we would obtain ^(A)A^ = with rrij the mass of this 

particle. However, for more complex X(x), such as the distance between two 
particles % and j, X(x) = jr, — r,-|, no simple analytic expression exists. The 



calculation of \0(X)X}^ can then be calculated by sampling a random set of 
A^mc velocities A as soon as a crossing is detected: 



e(x)x) = 





1 , „ Ef MC e(X)i 




1 

|A(xiAt;Ai)|_ 





(A.7) 



where i runs over all MD crossings with interface Ai, X(xi At ; Ai) is the MD 
crossing velocity through Ai, j runs over the N MC 'artificial' velocities A that 
are taken from a proper distribution P(X\x iAt ). For NVT simulations without 
additional constraints this distribution P(X\x iAt ) does not depend on the phase 
point x iAt and we can simply sample X({p}) where the momenta {p} defining A 
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are taken from a Gaussian distribution For NVE simulations, the distribution 
P(\\x iAt ) does depend x iAt and we have to change all momenta and distribute 
them on the hypersphere defined by the kinetic energy K = E — V(xiAt', Ai) 
with E the total energy and V(xiAt] Ai) the total potential energy at the 
crossing point. The proper sampling of momenta distributions in the presence 
of linear constraints, such as linear and angular momentum is explained in 

Ref. [37]. Clearly, if P(X\x iA t) = S(X - \(x iAt ; Ai)), Eq. (A. 7) would be equal 

to N+/J2i*\HxiAt, Ai)| _1 leaving Eq. (A.l) identical to Eq. (28) from which 
we started. 



B TIS shooting acceptance criterion for stochastic dynamics 



Although we assume throughout the paper that the equations of motion were 
deterministic, it is sometimes useful to implement some stochasticity into the 
dynamics, or consider completely stochastic equation of motion such as Brow- 
nian Dynamics [14,48,31]. Quantities like Jia( x o) are > then, no longer just 1 or 
but turn into probabilities with a fractional value. Moreover, for stochastic 
dynamics it is not trivial whether we are allowed to use the path that generated 
Xq as our instrument to search for the new phase point Xq . In this appendix 
we derive the acceptance probability for the shooting algorithm for arbitrary 
dynamics along the same lines as in Ref. [14]. At start, we try to be as general 
as possible making the least possible assumptions on the type of dynamics or 
on whether the system is in equilibrium or not. For this purpose, it is most con- 
venient to use the path space description, instead of phase space. The weight 
or probability density P[x] for a single path x = {x_ T t At , . . . ; x ; . . . , x +T f At } 
is than not only determined by the distribution p(x ) of x , but also by the 
probabilities of arriving along this precise route from x_ T b At in the past and 
continuing upto x +T j At in the future. 



P[x] = p{x ) n p{x(i + i) At <- X iAt ) nK^i-^At -> XiAt) (B.l) 
j=-l i=l 

where p(x — > y) is the forward transition probability (more accurate: proba- 
bility density) to go from x to y and p(y <— x) is the probability that, if the 
system is at y, it came from x in the past. Here, p(xo) is not necessarily the 
Boltzmann distribution or even a distribution in equilibrium. It is applicable 
to all systems that (at least to some approximation) are described by a steady 
state. We can express p(y <— x) in terms of forward transition probabilities as 



= p(x)p(x -> y) = p{x)p{x -> y) 
J (ix'p(x')p(x' -> y) p(y) 
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where / dx' p(x')p(x' — > y) — p(y) results from the steady state behavior. Using 
this relation, one can show that Eq. (B.l) is exactly equal to the probability 
of the first point p(x_ T t At ) times the forward evolution probabilities 

T /_! 

V[x] = p(x_ TbAt ) Y[ p(x iAt -> X( i+ i)At) (B.3) 

i=-r b 

which is identical to the weight for a path in the TPS-ensemble [14] with 
x -r b At instead of x . Note that so far, we have assumed nothing about the 
nature of dynamics (irreversible or reversible, stochastic or deterministic). 
When restricted to the TIS ensemble for interface i, the probability density of 
a path can be written as 

P Ai [x] = fe i (x)7'[x]/Z(A i ) (B.4) 

where hi is unity if the path goes from A , crosses \ and ends either at 
or goes back to A . Otherwise it is zero. The normalizing factor Z(X i ) equals 

Z{\) ee Jvxhi(x)V[x] (B.5) 

where the integral is taken over all possible paths x of all lengths, starting in 
all possible initial conditions X- Tb . Note that, contrary to TPS, Eq. (B.3) and 
Eq. (B.4) are not directly related to the relative probabilities of all paths in 
the TIS ensemble. This is a result of the path ensemble containing paths of 
different lengths. Eq. (B.3) turns into a true probability only when multiplied 
with the infinitesimal volume element in path space £>x = Hl = ^ T b dx iAt ~ 
dx N . Hence, a long path has an infinitely smaller probability than a shorter 
one for stochastic dynamics. Therefore, the concept of path space may sound 
peculiar for TIS. Still, it is instrumental to derive proper acceptance rules for 
TIS obeying detailed balance. 

When performing the random walk in the TIS path space using the shooting 
algorithm, the detailed balance condition is 



P gcn [x<°> -> x^] P acc [x W -> xW] 7\ [x W] 

P gen [xW - x(o)] P acc [xW - xW] V Xi [xW] 1 -° j 

where o and n, denote the old and new path respectively. The usual Metropolis 
acceptance rule is then 



mm 



P[X^] P g en[x( n ) -> Xg| 

' 7>[x(°)] P gen [x(°) -> xW] 



(B.7) 



Note that this rule only applies at the trajectory space level, it has nothing 
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to do with whether the underlying dynamics is stochastic or deterministic, or 
even reversible or irreversible. 

The generation probability to create a new path from an old path using the 
shooting move is given by 



P gcn [ X (°) - X W])] = ^^en[x (0) - xH]P g b cn [ X (°) - X<»>] (B.8) 

where is the chance to choose the shooting point r' with — r 6(o ) < r' < 

t^°^) at the old path, P(5p) the chance to select the randomized momenta 
displacements. As Sp is normally taken from a symmetric distribution, hence 
P(Sp) = P(—5p), this term will cancel in Eq. (B.7). The last two factors in 
Eq. (B.8) are the probabilities to generate trajectories from the shooting point 
point t' with the new momenta. These generation probabilities are given by 
the underlying dynamics used to generate the trajectories. If one starts from 
a shooting point at r' on the old existing path (with — rf^ < t' < rj°' ) ) the 
generation probability for the forward segment is 



-i 



^en[x (0) - XW] = J] P (*S "> 4l)At) ■ (B-9) 
i=r' 

This generation probability is exactly identical to the weight of the forward 
segment. The integration of the backward segment is not always trivial es- 
pecially when dealing with irreversible processes. However, in general, when 
reversible dynamics is applied, the backward segment is obtained by revers- 
ing the momenta, integrating forward in time and reversing the momenta 
again [14]. Accordingly, the backward segment's generation probability equals 

J£Jx (o) - xW] = if p (x$ 1)M - x\t) ■ (B.IO) 



where x = {r, —p} for a phase point x = {r, p}. Using these generation prob- 
abilities and the path weight Eq. (B.3) the factor within the min function of 
Eq (B.7) can be written as 



P[xW]P gcn [x(°)^x(°)] P [x_ r ]N^ 

P[x(°)]P gen [x(°) -> xW] p[x_ r(o) }N^ 1 ' 



r'-l r>[r (n) — ► T (n) 1 r'-l ^J ) _, ~(°) 

n P[ x iAt _ x (i+l)At\ -r-r FL x (i+l)At X »At 
r-(n) I(nh 11 ~T~(oj (S) 

=-r 6 ( n ) Pi X (i+l)At ~ ¥ x iAt\ i=-r 6 (o) Pl X iAt ~^ x (i+l)At 



where the generation probability and the weight of the forward parts of the 
trajectories have canceled each other. Eq. (B.3) simplifies tremendously for 
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dynamics that obey the microscopic reversibility condition [14] 

p(y^x) p(x)' 

This condition can be seen as a special case of Eq. (B.2) for time-reversible 
dynamics with p(y <— x) = p(y — > x), but is very general and valid for a broad 
class of dynamics applying to both equilibrium and non-equilibrium systems. 
As result, if the microscopic reversibility condition (B.12) is satisfied and thus 
p(y — > x) — p(y <— x), almost all terms in Eq. (B.ll) cancel, except for the 
steady state distributions p{x T i) of the shooting points. 

P acc [x(°) -^xW] =^(x( n ))min 

which is exactly the same as for deterministic dynamics. 

This is an important result as it allows to perform the acceptance/rejection 
rule (step 3 of Sec. 3.3) for the new momenta at the shooting point even if the 
energy along the path changes giving a different weight to p(aVAt) as to p(x ). 
The ratio /N^ in Eq. (B.13) can, of course, not be known in advance at 
the shooting point. However, this is effectively circumvented by defining N max 
at step 4 of the shooting algorithm in Sec. 3.3 leaving only the p{x^)) / p{x^)) 
term for the acceptance rule (step 3). 
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